//========================================================================================
// HashSamplingOverlapper.cc
//
// ----------------------------------------------------------------------------
// $Id: $
//------------------------------------------------------------------------------
// $Log: $
//
//
// KJD: Sometimes interval hits indicate span greater than source size. 
// Why does this happen, and what should be done about it?  For now, I truncate. 
//
//========================================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8.........9
//

#include "HashSamplingOverlapper.hh"

extern LogStream ls,ls2,ls3;
extern bool gDEBUG;

//---------------------------------------------------
// * ComputeGlobalRepeatStatistics
//---------------------------------------------------
// Pick some random samples and search them against the search engine
// to determine the distribution of hits for randomly chosen 
// samples.  Samples that hit significantly above the mean 
// can later be filtered as being repeats. 
//
//
void
HashSamplingOverlapper::ComputeGlobalRepeatStatistics(SequenceVector &Queries,
																		HashSampleSearchEngine &HSS)
{
  cout<<"\tCompute repeat statistics..."<<flush;

  RandomNumberGenerator RNG;
  SeqSampler SS;
  SequenceVector Samples;

  int numQueries = Queries.size();
  vector<int> HitsBySample;

  const int NumStatQueries = 2000; // Number of queries to use for repeat stats. 
  Samples.resize(NumStatQueries);
  HitsBySample.resize(NumStatQueries);
  
  //ls<<" NumStatQueries="<<NumStatQueries<<" SampleSize="<<mSampleSize<<flush;
  // Generate a bunch of random samples from the queries
  for(int i = 0;i< NumStatQueries;i++){
	 int qNum = (int)RNG.Uniform(0,numQueries-1);
	 // DEBUG  Look for zero length queries. 
	 for(int q = 0;q < Queries.size();q++){
		if (Queries[q].size() == 0){
		  cout<<"HEY.  Queries[q].size=0 q="<<q<<endl;
		  Queries[q-1].Print();
		  Queries[q].Print();
		  Queries[q+1].Print();
		  exit(1);
		}
	 }

	 // DEBUG  Testing for less than sample size since we're about to grab one.
	 if (Queries[qNum].size()<mSampleSize){
		cout<<"!!HEY!!! Queries[qNum].size() < mSampleSize"<<endl;
		cout<<"Queries[qNum].size()="<<Queries[qNum].size()<<endl;
		cout<<"mSampleSize="<<mSampleSize<<endl;
		cout<<"qNum="<<qNum<<" numQueries="<<numQueries<<endl;
		cout<<"Queries.size()="<<Queries.size()<<endl;
		cout<<"NumStatQueries="<<NumStatQueries<<endl;
		for(int q = 0;q < Queries.size();q++){
		  cout<<"Queries["<<q<<"].size="<<Queries[q].size()<<endl;
		}
		Queries[qNum-1].Print();
		Queries[qNum].Print();
		Queries[qNum+1].Print();
		exit(1);
	 }
	 SS.GetRandomSample(mSampleSize,Queries[qNum],Samples[i]);
	 if (Queries[qNum].size() == 0) cout<<"ERROR Query size = 0"<<endl;
	 if (Samples[i].size() == 0) cout<<"ERROR Samples[i] = 0"<<endl;

	 HitsBySample[i] = 0;
  }

  // Search them against the sources encoded in CST
  //SearchResultVec Results;
  SampleSearchResultVec Results;
  HSS.SearchWithSamples(Samples,Results);


  // Add up the number of hits for each sample
  long TotalHits = 0;
  for(int r = 0;r< Results.size();r++){
	 //	 if (Results[r].mSourceLoc >=0){  // KJD redundant
	 HitsBySample[Results[r].mSampleIdx]++;
	 TotalHits++;	
  }
  //ls<<"TotalHits "<<TotalHits<<endl;
  mMeanHitsPerSample = (float)TotalHits/(float)NumStatQueries;

  float totalDeviation = 0;
  for(int i = 0;i< NumStatQueries;i++){
	 double diff = HitsBySample[i] - mMeanHitsPerSample;
	 double deviation = diff*diff;
	 totalDeviation+=deviation;
  }
  double variance = (1.0/(NumStatQueries-1.0))*totalDeviation;
  mStDevHitsPerSample = sqrt(variance);
  // Actual cutoff for repeat samples.
  mRepeatThreshold = (int)((mStDevRepeatCutoff*mStDevHitsPerSample)+
									mMeanHitsPerSample); 
  
  cout<<"done."<<endl;

  cout<<"\tMeanHitsPerSample="<<mMeanHitsPerSample;
  cout<<" StDevHitsPerSample="<<mStDevHitsPerSample;
  cout<<" RepeatThreshold = "<<mRepeatThreshold;
  cout<<endl;
}

//---------------------------------------------------
// * ComputeSourceGroupSizes
//---------------------------------------------------
//
void
HashSamplingOverlapper::ComputeSourceGroupSizes()
{
  
  // figure out how much RAM we have, average size of a source, 
  // and come up with an average number of sources in each ST. 
  // This should accomodate worst cases. 
  float temp =((float)mpSources->size()/(float)mSourcesInGroup);

  int itemp = (int)temp;

  if (temp != itemp){
	 mNumSourceGroups = itemp+1;
  }else{
	 mNumSourceGroups = itemp;
  }
}

//---------------------------------------------------
// * AddToGraph
//---------------------------------------------------
// NOTE: OverlapGraph is assumed to already have space for 
// all of the nodes!
// 
void
HashSamplingOverlapper::AddToGraph(int QIdx,vector<float> &Overlaps,
											  AdjacencyListGraph & OverlapGraph)
{
  // Name the node with the query name
  OverlapGraph[QIdx].mName = (*mpQueries)[QIdx].mName;
  
  // For each overlap, if the overlap is large enough, put it in 
  // the OverlapGraph.  
  for(int i = 0;i< Overlaps.size();i++){
	 if (Overlaps[i] >= mMinOverlapForGraph){
		OverlapGraph[QIdx].AddEdge((*mpSources)[i].mName,Overlaps[i]);
	 }
  }
}

//---------------------------------------------------
// * ComputeHitsAndIntervalHitsOnSources
//---------------------------------------------------
// !! Must be called after ComputeSampleHitStats
//
void
HashSamplingOverlapper::ComputeHitsAndIntervalHitsOnSources(SampleSearchResultVec &SearchResults,
																				SampleHitStats &Stats)
{
  if (SearchResults.size() == 0) return;

  // Compute some statistics about hits for each sample
  for(int r = 0;r < SearchResults.size();r++){
	 Stats.mTotalHitsBySample[SearchResults[r].mSampleIdx]++;
  }

  SearchResults.SortBySourceLoc();
  SearchResults.SortBySourceIdx();
  
  int LastSampleIdx = SearchResults[0].mSampleIdx;
  int LastSourceLoc = SearchResults[0].mSourceLoc;

  // KJD r=0 correct start point for interval??  
  for(int r = 0;r< SearchResults.size();r++){		

	 // Search engine has been changed so that SourceIdx is relative to global source
	 // idx, not to the local results. 
	 int SourceIdx = (int) SearchResults[r].mSourceIdx;
	 int SampleIdx = (int) SearchResults[r].mSampleIdx;
	 int SourceLoc = (int) SearchResults[r].mSourceLoc;
		
	 // Samples which hit more than the derived repeat threshold are not counted
	 if (Stats.mTotalHitsBySample[SampleIdx] < mRepeatThreshold){
		(*Stats.mTotalHitsOnSource).IncrementOrAdd(SourceIdx,1); // Adds 1 if source in hash, creates it if not.
		(*Stats.mLastHitOnSource).AddOrChange(SourceIdx,SourceLoc);  // Adds if source not in hash, changes if it is. 
	 }else{
		Stats.mRepeatFilterCount++;
	 }
		
	 // The interval that these two samples are apart
	 long interval = SourceLoc - LastSourceLoc;
	 long targetInterval = (SampleIdx - LastSampleIdx)*mSampleSpacing;
	 
	 // Need something to check that samples are monotonic
	 
	 // If the intervals are close enough, record an interval hit
	 long intervalDiff = 	 abs(abs(interval)-abs(targetInterval));

	 if (intervalDiff < mSampleSize){
		(*Stats.mIntervalHitsOnSource).IncrementOrAdd(SourceIdx,1);
		(*Stats.mIntervalSpanOnSource).IncrementOrAdd(SourceIdx,abs(targetInterval));
	 }
	 LastSampleIdx = SampleIdx;
	 LastSourceLoc = SourceLoc;	
  }
}



//---------------------------------------------------
// * FilterInteriorSingletons
//---------------------------------------------------
// !! Must be called after ComputeHitsOnSource or ComputeHitsAndIntervalHits
// 
// Removes singleton hits that are not on the edges of reads. 
// 
void 
HashSamplingOverlapper::FilterInteriorSingletons(SampleHitStats &Stats)
{
  // KJD BUG:  Sometimes s is out of range of mpSources!!!

  // Filter interior singletons. Not efficient, but hey...
  for(int s=0;s<mpSources->size();s++){
	 // KJD fix so that it doesn't create an element. 
	 // If the hit is a singleton
	 if ((*Stats.mTotalHitsOnSource)[s] == 1){
		int idx;
		if ((*Stats.mLastHitOnSource).Contains(s)){
		  idx = (*Stats.mLastHitOnSource)[s];
		}else{
		  idx = -1;
		}
		// If the hit is in the interior
		if ((idx > 50) && (idx < (mpSources[s].size()-50))){
		  (*Stats.mTotalHitsOnSource)[s] = 0;
		} 
	 }
  }
}

long bugStartTime;
double bugCumulativeTime;
double bugElapsedTime;
double bugCount = 0;

//---------------------------------------------------
// * EstimateOverlaps
//---------------------------------------------------
//
void
HashSamplingOverlapper::EstimateOverlaps(vector<float> &Overlaps,
													  SampleHitStats &FwdStats,
													  SampleHitStats &RevStats)
{

  
  //for(register int s = mSourceStart;s<=mSourceEnd;s++){
  StatHash::iterator it;
  for(it = (*FwdStats.mIntervalSpanOnSource).begin();it != (*FwdStats.mIntervalSpanOnSource).end();it++){
	 bugCount++;

	 int revISpan;
	 int key = (*it).first;  // first == key (i.e. source Idx)
	 int fwdISpan = (*it).second; // Get forward interval span (second == value)
	 // Get corresponding reverse interval span, if there are any
	 if (!(*RevStats.mIntervalSpanOnSource).Lookup(key,revISpan)){
		revISpan = 0;
	 }

	 // Assign the overlap to whichever has the greater number of 
	 // interval spans, either the forward or the reverse. 
	 if (fwdISpan > revISpan){
		Overlaps[key] = fwdISpan;
	 }else{
		Overlaps[key] = revISpan;
	 }	 
  }

  //Overlaps[s+mSourceStart] = mIntervalHitsOnSource[s]*mSampleSpacing;
  //Overlaps[s+mSourceStart] = mTotalHitsOnSource[s]*mSampleSpacing;
  //Overlaps[s+mSourceStart] = mTotalHitsOnSource[s];
}

//---------------------------------------------------
// * ConvertOverlapsToFraction
//---------------------------------------------------
//
void
HashSamplingOverlapper::ConvertOverlapsToFraction(vector<float> &Overlaps,
																  SequenceVector &Seqs)
{
  for(int s = mSourceStart;s<=mSourceEnd;s++){
	 double SourceSize = Seqs[s].size();
	 double fraction = Overlaps[s]/SourceSize;
	 int fractionI = int(fraction*100);
	 if (fractionI > 100){
		fractionI = 100;
		//cout<<"ERROR"<<endl;
		//cout<<"s="<<s<<"\t"<<"Over="<<Overlaps[s]<<"\tSource="<<SourceSize<<endl;
	 }

	 Overlaps[s] = fractionI;
  }
}

void
HashSamplingOverlapper::ComputeOverlapsWithSources(SampleSearchResultVec &SearchResults,
																	SampleHitStats &Stats,
																	SequenceVector &Seqs,
																	vector<float> &Overlaps)
{
  if (SearchResults.size() == 0) return;

  // Compute some statistics about hits for each sample
  for(int r = 0;r < SearchResults.size();r++){
	 Stats.mTotalHitsBySample[SearchResults[r].mSampleIdx]++;
  }

  SearchResults.SortBySourceLoc();
  SearchResults.SortBySourceIdx();
  
  int LastSampleIdx = SearchResults[0].mSampleIdx;
  int LastSourceLoc = SearchResults[0].mSourceLoc;

  long IntervalHitsOnSource = 0;
  long IntervalSpanOnSource = 0;

  // KJD r=0 correct start point for interval??  
  //for(int r = 0;r< SearchResults.size();r++){		
  int r = 0;
  do{
	 // Search engine has been changed so that SourceIdx is relative to global source
	 // idx, not to the local results. 
	 int CurrentSourceIdx = (int) SearchResults[r].mSourceIdx;
	 // As long as we are on the same source, accumulate hits and interval hits. 
	 do{
		int SampleIdx = (int) SearchResults[r].mSampleIdx;
		int SourceLoc = (int) SearchResults[r].mSourceLoc;
		
		// Samples which hit more than the derived repeat threshold are not counted
		if (Stats.mTotalHitsBySample[SampleIdx] < mRepeatThreshold){
		  long interval = SourceLoc - LastSourceLoc;
		  long targetInterval = (SampleIdx-LastSampleIdx)*mSampleSpacing;
		  long intervalDiff = abs(abs(interval)-abs(targetInterval));
		  
		  if (intervalDiff < mSampleSize){ // KJD arbitrary heuristic
			 IntervalHitsOnSource++;
			 IntervalSpanOnSource+=targetInterval;
		  }
		  LastSampleIdx = SampleIdx;
		  LastSourceLoc = SourceLoc;
		}
		r++;
		if (r >= SearchResults.size()) break; // End this loop if we have exhausted results
	 }while(SearchResults[r].mSourceIdx == CurrentSourceIdx);
	 // Now that we've looked at all the samples for the given source, estimate 
	 // the overlap between that source and the query and save it in Overlaps
	 // For now, the estimate is simply the IntervalSpan
	 double OverlapFraction = (double)IntervalSpanOnSource/(double)Seqs[CurrentSourceIdx].size();
	 int fractionI = int(OverlapFraction*100);
	 if (fractionI > 100) fractionI = 100;
	 // If called with both fwd and rev results, want to save the larger of the two
	 if (Overlaps[CurrentSourceIdx] < fractionI) Overlaps[CurrentSourceIdx] = fractionI;

	 IntervalHitsOnSource = 0;
	 IntervalSpanOnSource=0;
  }while(r < SearchResults.size());

}


//---------------------------------------------------
// * ComputeOverlapsForQuery
//---------------------------------------------------
// Computes the overlaps for a single query against all of the
// sequences in the CST.   
//
// KJD Note: Queries are assumed to be the same as sources here
// because Queries is passed into ComputeOverlapsWithSources as
// a proxy fo Sources KJD Be sure to change if source != query. 
//
void 
HashSamplingOverlapper::ComputeOverlapsForQuery(int QueryIdx,
													 HashSampleSearchEngine &HSSE,
													 SequenceVector &Queries,
													 AdjacencyListGraph &OverlapGraph)
{
  vector<float> Overlaps(mpSources->size(),0.0);

  // Clear the stats from the previous read. 
  //FwdStats.Clear();
  //RevStats.Clear();

  // Compare query samples to current sources
  SampleSearchResultVec FwdSearchResults,RevSearchResults;
  StartTiming(bugStartTime);
  long FwdNumSamples = HSSE.Search(Queries[QueryIdx],FwdSearchResults);
  long RevNumSamples = HSSE.SearchRev(Queries[QueryIdx],RevSearchResults);
  bugElapsedTime = EndTiming(bugStartTime);
  bugCumulativeTime+=bugElapsedTime;

  // KJD: Shouldn't do this before call ComputeOverlapsForQuery?
  // Return if sequence was too small to sample. 
  if ((FwdNumSamples <= 0) || (RevNumSamples <=0)) return;

  // Look at FORWARD samples ------------------
  FwdStats.InitSampleStats(FwdNumSamples); 

  // Sometimes search returns no results, so there is nothing to do. 
  if (FwdSearchResults.size()>0){
	 //ComputeHitsAndIntervalHitsOnSources(FwdSearchResults,FwdStats);

	 ComputeOverlapsWithSources(FwdSearchResults,FwdStats,Queries,Overlaps);
	 // if (mbFilterInteriorSingletons) FilterInteriorSingletons(FwdStats); 
  }

  // Look at REVERSE samples ------------------
  RevStats.InitSampleStats(RevNumSamples);

  // Sometimes search returns no results, so there is nothing to do. 
  if (RevSearchResults.size()>0){	
	 //ComputeHitsAndIntervalHitsOnSources(RevSearchResults,RevStats);
	 ComputeOverlapsWithSources(RevSearchResults,RevStats,Queries,Overlaps);
	 //if (mbFilterInteriorSingletons) FilterInteriorSingletons(RevStats);
  }
  
  // Use those statistics to estimate overlap amounts, choosing best overlap from 
  // Fwd or Rev samples.  
  
  //EstimateOverlaps(Overlaps,FwdStats,RevStats);

  // Overlaps are initialy in bp.  
  //ConvertOverlapsToFraction(Overlaps,*mpSources);

  // Add this group of overlaps to the growing OverlapGraph
  AddToGraph(QueryIdx,Overlaps,OverlapGraph);
}


//---------------------------------------------------
// * ComputeOverlapGraph
//---------------------------------------------------
// Computes the complete (sparse) overlap graph for the given 
// set of sources and queries. 
//
//
void
HashSamplingOverlapper::ComputeOverlapGraph(AdjacencyListGraph &OverlapGraph)
{

  bool bFirstTime = true;
  FwdStats.mRepeatFilterCount = 0;
  RevStats.mRepeatFilterCount = 0;
	 
  // Initialize stats vectors. 
  //InitSourceStatVectors(mpSources->size(),FwdStats); 
  //InitSourceStatVectors(mpSources->size(),RevStats); 

  // Add nodes for all of the queries 
  OverlapGraph.resize(mpQueries->size());

  // For each source group (i.e. suffix tree or hash of reads).
  ComputeSourceGroupSizes();
  mSampleSpacing = (int)(1000.0/(float)mSamplesPerKbp);
  for(int SourceGroup = 0;SourceGroup < mNumSourceGroups;SourceGroup++){

	 bugStartTime = 0;
	 bugCumulativeTime = 0;
	 bugElapsedTime = 0;

	 cout<<"SourceGroup#: "<<SourceGroup<<" "<<flush;

	 // Create a search engine for the sources in this group. 
	 HashSampleSearchEngine *pHSSE = new HashSampleSearchEngine();
	 pHSSE->mSampleSize = mSampleSize;
	 pHSSE->mSampleDensity = mSamplesPerKbp;
	 pHSSE->mTotalSamples = 0;
	 pHSSE->mMaxLocsPerMer = mMaxLocsPerMer;
	 
	 mSourceStart = (SourceGroup*mSourcesInGroup);
	 mSourceEnd = (SourceGroup+1)*mSourcesInGroup-1;
	 if (mSourceEnd >= mpSources->size()) mSourceEnd = mpSources->size()-1;
	 int NumSources = mSourceEnd - mSourceStart + 1;

	 cout<<"  (SourcesInGroup: "<<NumSources;
	 cout<<"  Source Start: "<<mSourceStart<<" SourceEnd: "<<mSourceEnd<<")"<<endl;

	 long StartTime;
	 StartTiming(StartTime);	 
	 cout<<"\tCreate source hash.."<<flush;
	 pHSSE->Create(*mpSources,mSourceStart,mSourceEnd);
	 float elapsedTime = EndTiming(StartTime);
	 cout<<"done. "<<elapsedTime<<" sec."<<endl<<endl;

	 // Compute some statistics about the distribution of repeat 
	 // samples (for filtering).
	 if (bFirstTime){
		StartTiming(StartTime);
		//cout<<"mpQueries->size()="<<mpQueries->size()<<endl;
		ComputeGlobalRepeatStatistics(*mpQueries,*pHSSE);
		float ElapsedTime = EndTiming(StartTime);
		cout<<"\t"<<ElapsedTime<<" seconds."<<endl<<endl;;
	 }

	 StartTiming(StartTime);
	 cout<<"\tCompute overlaps "<<mpQueries->size();
	 cout<<" queries vs. "<<NumSources<<" sources.."<<flush;
	 // Compare every query to the sources in HSSE
	 for(int QueryIdx = 0;QueryIdx < (*mpQueries).size();QueryIdx++){
		if (!(QueryIdx%1000)) cout<<"."<<flush;
		ComputeOverlapsForQuery(QueryIdx,*pHSSE,(*mpQueries),OverlapGraph);
	 }
	 float ElapsedTime = EndTiming(StartTime);
	 cout<<"done. "<<endl;
	 cout<<"\tQuery Samples: "<<pHSSE->mTotalSamples;
	 cout<<"   TotalTime: "<<ElapsedTime<<" seconds.";
	 cout<<"   SearchTime: "<<bugCumulativeTime<<" seconds."<<endl;	 
	 cout<<endl;
	 cout<<"\t MaxLocPerMer="<<mMaxLocsPerMer<<"  ";
	 cout<<" MaxLocCount="<<pHSSE->mbugMaxLocsCount<<endl;
	 delete pHSSE;
  }
}

